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Abstract Modelling of turbine blade-induced turbulence (BIT) is discussed within 
the framework of three-dimensional Reynolds-averaged Navier-Stokes (RANS) ac- 
tuator disk computations. We first propose a generic (baseline) BIT model, which 
is applied only to the actuator disk surface, does not include any model coefficients 
(other than those used in the original RANS turbulence model) and is expected to 
be valid in the limiting case where BIT is fully isotropic and in energy equilibrium. 
The baseline model is then combined with correction functions applied to the region 
behind the disk to account for the effect of rotor tip vortices causing a mismatch of 
Reynolds shear stress between short- and long-time averaged flow fields. Results are 
compared with wake measurements of a two-bladed wind turbine model of Medici 
and Alfredsson [Wind Energy, Vol. 9, 2006, pp. 219-236] to demonstrate the capa- 
bility of the new model. 



1 Introduction 

Modelling of turbulent mixing behind turbines is an important area of research for 
the wind and tidal power industries. Whilst recent eddy-resolving computations 
(such as large-eddy simulations and detached-eddy simulations) using a so-called 
actuator line method have given valuable insight into time-dependent features of 
wind/tidal turbine wakes (e.g., 12), it is still of great importance to develop a low- 
order model that satisfactorily predicts the characteristics of time-averaged flow 
around turbines at much lower computational cost. Such a low-order model would 
be useful for the design of wind and tidal power farms for the future. 
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The objective of the present study is to develop a reasonably generic low-order 
model that describes time-averaged flow around a horizontal-axis wind/tidal turbine 
of various designs and operating conditions. Our current model is based on the three- 
dimensional Reynolds-averaged Navier Stokes (RANS) equations coupled with an 
actuator disk model. To account for the effect of turbulence (or velocity fluctuations) 
induced by turbine blades — we refer to this as blade-induced turbulence (BIT) in 
this paper — we first introduce a simple baseline model, which requires the specifi- 
cation of BIT energy and scale but does not require any additional model coefficients 
(other than those used in the original RANS turbulence model). The baseline model 
is then combined with correction functions to account for the effect of rotor tip vor- 
tices. Details of the model are described below, followed by the results of a model 
validation study. 



2 Baseline Model 

The baseline model used in this study is the one recently proposed by the authors 
Q. The governing equations of the flow are the three-dimensional incompressible 
RANS equations, where the Reynolds stress tensor, —u\u'-, is modelled using the 
standard k-e model of Launder and Spalding |2|. 

A simple actuator disk concept is used, where a turbine rotor is modelled as a 
stationary permeable disk of zero thickness placed perpendicular to the incoming 
flow. The effect of the rotor on the mean flow is considered as a loss of streamwise 
(or axial) momentum at the disk plane. The change in momentum flux (per unit 
disk-area and per unit fluid-density) is locally calculated as Sy = \KUj, where t/j 
is the local streamwise velocity at the disk plane and K is a. momentum loss factor 
(constant over the disk) to determine the thrust acting on the disk. 

It is known that this actuator disk model (without taking account of BIT) yields 
much weaker turbulent mixing behind the disk compared to that measured behind a 
turbine rotor In the baseline BIT model we assume that, at the disk plane, virtual tur- 
bine blades generate turbulence that is characterised by its turbulent kinetic energy 
kj, and dissipation rate £/,. Two physical parameters are introduced to determine the 
values of kf, and e^: (i) the ratio of the energy converted to BIT to that removed from 
the mean flow at the disk plane, j3, and (ii) a representative length scale for BIT, //,. 
Note that j3 and li, are model variables rather than model coefficients, as they de- 
pend on actual turbine design and operating conditions, ki, is then calculated (locally 
over the disk plane) as = fiSu = jj3KUj, whereas ei, is estimated (based on the 

3 3 

high Reynolds number equilibrium hypothesis) as = C^k^ /4, where = 0.09 
(following the standard k-e model). 

To account for the combined effect of BIT and the turbulence coming from up- 
stream of the turbine, we further assume that: (i) the disk plane is a special internal 
boundary that does not allow the transport of A: or e via their diffusion (so that the 
transport of k and e through the disk plane is only via the convection from upstream 
to downstream) and (ii) at the disk plane, the BIT of kj, and £/, is mixed (in the time- 
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averaged sense) with the ambient (or upstream) turbulence of ka and e„, resulting in 
the mixed turbulence of fe,„ and e,„ just downstream of the disk. Here k,„ is calcu- 
lated as k,n = ka +k[„ whereas e„, is estimated by assuming the conservation of the 
time-integral of linearly decaying turbulent kinetic energy, i.e., fe„,T,„ — kaTa + k^Tb, 
where Tq = ka/£a, ib = kb/Eb and Tm = km/e,,, represent the initial eddy turnover 
time or lifetime for the ambient, blade-induced and mixed turbulence, respectively 
||5l. Eventually, the changes in k and e to be added to their transport equations at the 
disk plane (per unit disk-area) are calculated as 

Sk = Ud {k,„ -ka) = Udkb , ( 1 ) 

{ka+kbf 



Se = Ud{£m — £«) = Ud 



(2) 



_{kl/ea) + {kl/eb) 

Generally, the two model variables, j3 and lb, should depend on the actual turbine 
design and may be given as functions of the distance, r, from the disk axis. In the 
present study, however, uniform values of )3 and lb are given either to the entire disk 
surface or to the disk edge region defined by {0.5d — Wedge) < < 0.5d, where d 
is the disk diameter and Wedge (= O.lii in this study) is the width of the disk edge 
region. For the latter case, )3 = is given to the rest of the disk surface. 



3 Tip Vortex Correction 

As will be shown later, the baseline model tends to yield too strong/fast turbulent 
mixing in the core region (r/d < 0.4) in the near wake, even when BIT is given only 
around the disk edge. This suggests that further modifications are required on the 
modelling of turbulent mixing in the near wake region. 

An important issue to be considered here is how large the Reynolds shear stress 
is (relative to the turbulent kinetic energy) behind a turbine, especially in the region 
where strong rotor tip vortices exist. It is known that for many two-dimensional 
turbulent shear flows where the production and dissipation of k are close to equi- 
librium, the value of u[u'^,/k is around 0.3 (when the mean shear dU/dy < 0) or 
—0.3 (when dU/dy > 0) and this is the basis of the turbulent viscosity constant 

= (0.3)^ = 0.09 used in the standard k-e model (see, e.g.. Pope 161). In the mean 
shear layer behind a turbine rotor, however, experimental results have shown that 
the magnitude of m^m' /A: can be significantly smaller than 0.3 [4\. This is most likely 
due to the effect of tip vortices causing a mismatch between short-time (or phase) 
averaged and long-time averaged wake shear profiles. Specifically, tip vortices peri- 
odically create "adverse" velocity gradient regions, d{u)/dr <Q (where (0) denotes 
a short-time average of over a time scale relevant to the passing-through time of 
each tip vortex), and hence "adverse" Reynolds shear stress regions, (m' m'^) > 0, 
whilst the long-time averaged wake shear direction is dU/dr > and therefore 
u'^u'^ < 0. This explains why the long-time-averaged wake shear behind a rotor does 
not diffuse quickly despite its high level of turbulent kinetic energy. 
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To account for the above effect of tip vortices within the framework of RANS 
actuator disk computations, we consider applying empirical correction functions to 
the so-called tip vortex region. The streamwise extent of the region will be given as 
a model variable, /tip, which may be further modelled as a function of turbine design 
and operating conditions in future studies. For the present study, we introduce three 
correction functions, /(- and f^i, into the standard k-e model: 

A ^^/^_£ (3) 

Dt dxj \ (Ji dxj J ' dxj ' 
De d fv, de\ —dUi . ^ e ^ ... 

- u'-u . ^r— /eiQi 7 - Q2 — , (4) 



Dt dxj \ <j£ dxj J ' dxj k k 

where = 0.09, Qi = 1 .44, Q2 = 1 -92, cta: = 1 .0 and CTe 1 .3 (following the stan- 
dard k-z model). To reduce the eddy viscosity and thus mitigate turbulent mixing in 
the tip vortex region behind the disk (located at x = 0), we model /|, as follows: 

fyi = \ for X < , /tip < , (6) 
= !-/,-( l-/^™i„) for 0<x</*p, (7) 

/^ = l-/,(l-/^,nin):r 1+cos- for Z*p<x</„p, (8) 

^ \ 'tip ^ 'tip / 

where /^min = 0.1 is given in this study, following experimental observations ||4l 
(note, however, that this is also a model variable and should generally depend on 
the rotor tip-speed ratio and the number of blades). /*jp is the streamwise extent of 
the region where does not change (corresponding to the region where tip vortices 
are stable); we assume Z*p = jZnp in the present model, fy is a damping (sinusoidal) 
function in the radial (r) direction, normalised such that = 1 at r = O.Sc/ (where 
the rotor edge is located) and monotonically decreases to = at r = and d. 

A major difficulty in this correction is that the change in v, due to the introduc- 
tion of affects not only the strength of turbulent diffusion but also the production 
of k and e. This necessitates the introduction of fj^ and f^i in Eqs (3) and (4), re- 
spectively; however, we assume /gi = 1 in this study for simplicity}^ Preliminary 

' Here the introduction of is justified since, in the tip vortex region, the production of k should be 
hnked to the mean of the magnitude of short-time averaged velocity gradient and the corresponding 
Reynolds stress rather than those for the long-time averaged flow field, i.e., /<- accounts for the 

difference between — kJk^.^ (known in the model) and | — (ujuy)^^! (unknown in the model). 
Meanwhile, physical interpretation of seems less straightforward since the production term 
of e includes e and k, which should also be divided into two components corresponding to the 
short-time and long-time averaged flow fields, respectively. This seems far beyond the capability 
of single-scale eddy viscosity models (perhaps we require a proper multi-scale model for such 
discussion; see, e.g., Sagaut et al. 1 7 1) and hence we assume /^i = 1 in this study. 
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computations with fk = I yielded too small k in the near wake, whereas those with 
fk — ^ /fji (i.e., fully cancelling out the effect of on the production of k) resulted 
in too large k in the near wake. To maintain the level of turbulent kinetic energy 
in the near wake comparable to that measured in the experiments |4 | (and also to 
minimise its sensitivity to the value of /up), we model fk as /t = 1 + 7a- ( 1 — ) /fji 
with the model coefficient 7. = 0.85 in this study (note that Yk=0 and 1 correspond 
to the two extreme cases, fk—^ and 1 / respectively). 



4 Model Validation 

A model validation study was performed for wind tunnel tests of a two-bladed wind 
turbine model of Medici and Alfredsson |l3]|4l. The tunnel test section is 1.2m wide 
and 0.8m high. The tunnel ceiling is adjustable, so that the freestream velocity, Ui^, 
is constant throughout the test section (when the turbine is not installed). The rotor 
diameter d is 0. 18m and its hub height from the floor is 0.24m. The hub diameter is 
20mm. Table 1 summarises experimental conditions; for this particular set of exper- 
iments, a turbulence-generating grid was installed, providing freestream turbulence 
(FST) intensities of 4.5% at x/d = (turbine location) and 2.5% at x/d = 9. Here 
the thrust coefficient Ct and the tip-speed ratio (TSR) are based on the freestream 
velocity and therefore slightly different from those reported earlier [3 1. 

For the computations, the turbine rotor was modelled as an actuator disk of zero 
thickness with a small non-permeable disk embedded at the centre to account for 
the effect of the rotor hub. The supporting tower of the turbine was not modelled. 



Table 1 Experimental conditions 



t/in 


TSR 


Ct 


Yaw angle 


FST intensity 


8.14 [m/s] 


3.87 


0.899 


[deg.] 


4.5% (x/d = 0), 2.5% (x/d = 9) 



Fig. 1 Contours of turbulent 
kinetic energy across a hor- 
izontal plane at hub height: 
(a) without BIT, (b) baseline 
model (for the entire actuator 
disk), (c) baseline model (only 
for 0.4 < r/d < 0.5), and (d- 
f) baseline model (only for 
0.4 < r/d < 0.5) with the tip 
vortex correction (hip/d = 1, 
3, 5 for d, e, f, respectively). 
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Fig. 2 Centreline stream wise velocity (left), streamwise velocity profiles at x/t/ = 3 (centre) and 
turbulent kinetic energy profiles at x/d = 3 (right): (a) without BIT, (b) baseline model (for the 
entire actuator disk), (c) baseline model (only for 0.4 < r/d < 0.5), (d) baseline model (only for 
0.4 < r/d < 0.5) with the tip vortex correction {kip/d = 3), and experimental data |3|. 



Computations were performed using FLUENT 12 (together with its User-Defined 
Function module). The solver is based on a finite volume method and is nominally 
2nd-order accurate in space. All computations were performed as steady state. 

Figures 1 and 2 present results for Ct = 0.899, j3 = 0.30 and 4 = O.ld (except 
for the case without BIT). The case without BIT yields too small turbulent kinetic 
energy and thus too slow wake recovery. The baseline model applied to the entire 
disk yields too large turbulent kinetic energy in the core region (r/d < 0.4) and 
hence too early wake recovery behind the disk. The baseline model applied only to 
the edge region still yields too strong turbulent mixing in the core region in the near 
wake, and the tip vortex correction provides better results in both near and far wake. 
Further validation studies are needed to examine the robustness of the model. 
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